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Abstract 

We  present  a  class  of  nonlinear  dynamic  viscoelastic  models  for  materials  subjected 
to  shear  stress.  The  model  equations  are  based  on  a  continuum  variation  of  a  reptation 
model  in  which  chemically  cross-linked  (CC)  systems  of  molecules  act  as  constraint 
boxes  per  unit  volume  for  physically  constrained  (PC)  systems  of  molecules.  Results 
from  validating  the  model  with  dynamic  shear  experiments  are  given  and  a  stability 
analysis  for  the  corresponding  linearized  systems  is  discussed. 


1  Introduction 

In  this  paper  we  derive  for  the  first  time  a  continuum  model  combined  with  molecular  based 
internal  dynamics  for  shear  deformations  in  viscoelastic  materials.  The  resulting  partial 
differential  equation  model  is  coupled  to  ordinary  differential  equations  for  internal  strains 
via  a  nonlinear  stick-slip  molecular  theory.  An  initial  step  toward  validating  the  models 
with  carefully  designed  experiments  is  discussed. 

Various  molecular  and  phenomenological  models  have  been  proposed  to  model  rubber  de¬ 
formations.  A  more  scientifically  satisfactory  model  can  be  obtained  by  combining  molecular 
and  continuum  mechanics  approaches.  The  fundamental  ideas  originate  in  the  presentations 
of  Doi  and  Edwards  [8]  and  Johnson  and  Stacer  [11]  as  extended  in  [4,  5].  To  make  the 
presentation  in  this  paper  self  contained  we  first  briefly  summarize  some  of  the  background 
material  presented  in  [3,  4,  5]. 

A  dynamic  model  for  rubber  viscoelasticity  has  been  proposed  in  [11]  based  on  a  con¬ 
tinuum  simulation  of  the  reptation  model  considered  in  [8].  In  step-strain  relaxations  of 
polymers  with  constraint  (stick-slip)  theory,  PC-molecules  deform  with  CC-molecules,  con¬ 
tract  and  creep  to  return  to  a  lower  energy,  higher  entropy  state.  As  a  result,  the  total 
energy  density  at  a  constant  strain  dissipates  in  time  and  viscoelastic  behavior  results. 

Such  a  linear  continuum  model  of  rubber  viscoelasticity  is  developed  in  [11]  by  considering 
the  CC-linked  molecules  as  providing  cells  or  boxes  with  entrapped  molecular  segments  and 
placing  a  unit  cell  or  box  at  each  point  of  the  rubber  continuum.  This  approach  was  extended 
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to  nonlinear  models  for  rubber  undergoing  tensile  deformations  in  [4,  5].  Here  we  consider 
the  case  where  the  CC-box  experiences  a  shear  stress  as  a  result  of  which  the  PC-molecules 
undergo  a  parallel  deformation.  In  our  formulation  the  PC-molecules  act  as  internal  variables 
in  the  relaxation  of  the  CC-box. 


2  Continuum  model  of  Johnson  and  Stacer 

When  the  Johnson-Stacer  [11]  CC-system  experiences  shear  deformations,  the  PC-system 
undergoes  deformation  along  with  it.  The  relaxation  of  the  PC-system  is  determined  by 
the  history  of  the  CC-system  along  with  thermodynamic  considerations.  On  the  other  hand 
the  PC-system  deformation  acts  as  an  internal  variable  affecting  the  relaxation  of  the  CC- 
system.  The  tendency  for  the  PC-moleculc  to  return  to  a  lower  energy  and  higher  entropy 
configuration  forms  the  basis  of  the  viscoelastic  model  in  the  deformation  process.  A  PC- 
moleculc  entrapped  in  a  CC-molecular  box  undergoes  a  deformation  in  the  direction  of  the 
shear  when  the  CC-box  experiences  shear.  After  the  initial  induced  stress,  the  PC-molecule 
relaxes  even  if  the  CC-box  remains  deformed.  Let  the  CC-constraint  tube  have  length  L(t) 
and  the  representative  entrapped  PC-moleculc  have  length  £(£).  The  model  in  [11]  assumes 
that  the  length  of  the  CC-tube  is  approximated  as  a  step  function  of  time  and  £(0)  =  L( 0). 
Suppose  instantaneous  step-strain  deformation  of  the  CC-system  results  in  the  PC-molecule 
having  length  £*.  A  model  of  the  PC-molecule  returning  to  its  original  contour  length  is 
given  by 

e(t)  =  £(0)  +  [t  -  ^(0)]  e~t/T 

If  the  initial  stretch  for  the  CC-box  is  A Ln,  relating  the  stretches  in  the  form  —  =  — 

AL0  L0 

leads  to  the  relation 

A<°  =  (|) 

Then  in  the  time  interval  to  <  t  <  t\  we  have 

Kt)  =  4  +  y-AL0e_^_to^r, 

where  r  is  the  relaxation  time  for  the  slip  motion  of  the  PC-molecule. 

Setting  t  —  ti,  we  have 


£(ti)  =  4  +  -T-A  L0e"(tl"to)/r, 
Aj 

Li  =  L0  +  AL0- 


Next  we  suppose  that  the  CC-box  is  subjected  to  an  instantaneous  shear  step  deformation 
A  Li  at  t  =  t\.  For  t  in  the  time  interval  (ti,t2)  we  have 

£(t)  =  4  +  [(^i  +  AA)-4]e-(t-tl)/T 

=  io  +  ^AL0e'(f-to)/T  +  -^A  Lie-{t~tl)/T . 

L  o  L  i 


2 


Similarly,  after  another  instantaneous  shear  step  deformation  AL2  at  t  —  £2,  for  £2  <  t  <  ts 
we  have 


£{t)  =  £0  +  ^AL0e_(t'fo)/r  +  ^-ALle~{-t~tl)/T 
Lq  L 1 

+  ^-A  L2e~{t-t2)/T. 

L‘2 

Following  [11]  and  taking  the  limit,  as  At*  — *  0,  of  a  succession  of  such  step  strains,  we  are 
led  to  the  formula 

l(t)=  1 »+  1‘tT)  (2.1) 

Jo  Lis)  ds 

for  the  length  of  the  PC-molecule. 

We  make  the  following  observation  for  the  linear  theory  that  will  be  subsequently  modified 
for  a  nonlinear  theory.  If  A  ^  and  A  Li  represent  instantaneous  step-stretches  at  t,  —  tt,  we 
have  assumed  the  fundamental  relationship 


A  Lt 


(2.2) 


where  ^  =  £(ti),  L *  =  L(ti). 


3  Shear  deformation  in  the  CC-box  and  deformation 
of  the  PC-molecule 

A  deformation  that  is  close  to  the  axial  extension  and  to  which  we  can  readily  apply  the 
continuum  ideas  of  Johnson  and  Stacer  is  generalized  simple  shear  as  depicted  in  Figure  1. 
We  suppose  that  the  unit  CC-box  undergoes  a  deformation  of  the  type  x  — ■>  x  +  u(y),  y  y, 

z  — »  z. 


y 


X 


> 


Figure  1:  Simple  shear. 
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The  configuration  gradient  [1,  2]  for  this  deformation  is  given  by 

(  1  u'(y)  0  \ 

A=  0  1  0 

\0  0  1J 

The  configuration  gradient  A  can  be  written  in  a  unique  way  as  a  product  of  a  stretch  tensor 
E  and  rotation  R  as 

A  =  ER, 

where  E2  =  A* A.  The  eigenvalues  £1,^2,  £3  of  A* A  are  given  by 

6  =  1  +  \  KM)2  +  j  [“'(</)]2 

=  1  +  <•'(</)  +  \  K(i/)]2  + 1  K(s/)]3  -  ■  ■  ■  (3.1) 

&  =  l  +  \{u'(y)?-u'(s,)sJ\  +  \[u'(s,tf 

=  1  -  <«'(</)  +  \  Kte)]2  - 1  K(!/)]3  +  ■  ■  ■ 

(3  =  1. 

We  note  that  £i£2  =  1  and  thus  also  ^^2^3  =  1-  For  small  deformations  we  set 
Aic  =  v7^  =  1  +  \u\v)  +  |  [u\y)f  +  ^  [■ u\y)f  T  •  •  • 

A2c  =  v7^  =  1  -  ^u\y)  +  ^  [u'(y)f  -  ^  [u'(y)f  +  . . .  (3.2) 

A3C  =  1- 


The  quantities  Aic,  A2c,  A3C  are  the  principal  stretches  for  the  CC-shear  deformations. 

Although  the  CC-box  experiences  shear  deformation,  the  entrapped  PC-molecules  expe¬ 
rience  deformations  that  can  be  simple  shear,  simple  elongation,  or  a  combination  thereof, 
depending  on  the  orientation  of  the  PC  molecule  in  the  entanglement  with  the  CC-box.  In 
any  case,  we  assume  that  the  PC  deformation  5up  is  in  the  same  direction  as  the  CC-box 
shear,  and  thus  the  principal  stretches  for  the  PC-molecules  may  be  assumed  to  obey  the 
relationships 

A  2P  =  A3P  =  _ —  (3.3) 

V  Aip 

To  develop  relationships  between  the  principal  stretches  A *c  and  A jP  we  use  (2.1)  and  (2.2). 
The  equation  (2.1)  suggests  that  we  write 


Aip  —  1  + 


Aip(s)  dXic(s)  _ 
Aic(s)  ds 


(t ~S)/Tds . 


(3.4) 


The  equation  (2.2)  further  suggests  that  we  assume  the  relationship 


OXjp 

d\jc 


\ 

A jc 


(3.5) 
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where  &l3  is  the  Kronecker  delta. 

In  (3.2)  we  assume  that  u  also  depends  on  time.  Thus, 

Aic  =  1  +  ^dyu(t,  y)  +  ^  [dyu(t,  y )]2  +  ^  [dyu(t,  y)]3  =F  •  •  • 

A2c  =  1  -  ^dyu(t,  y)  +  ^  [dyu(t,  y )]2  -  [dyu(t,  y)f  ±  . . .  (3.6) 

A3C  =  1- 


Using  (3.6),  we  obtain  from  (3.4)  the  equation 

Aip  =  1  +  ^  Aip(s)  ^1  -  yu{s,  y)  +  [dyu(s,  y)f  +  . .  • 

11  \ 

-d,dyu(s,y)  +  -dyu  (s,y)d,d,u(s,y)  +  ...]•  e_,l_*)/Tds. 


(3.7) 


From  (3.7)  we  hnd 

Xlp  ~  1  +  L  XlP ^  l1  _  \dyU^S,y^)  '  \^dyU(yS}y^e~{t~S)/Tds'  (3-8) 

Setting  Aip  ~  1  +  8up  where  5up  =  Sup(t,y),  we  see  from  (3.8)  that 

1  1  ft 

Sup(t,y)  «  -dyu(t,y)-—  dyu(s,y)e^{t~s)/rds 

2  2r  J 0 

^dsd yu(s,  y)e~(t~s),Tds  -  ^dyu( 0,  y)e~t/r.  (3.9) 

Next,  let  (assume  that  dyu(0,y)  =  0) 

1 

w  —  8up -dyu(t,  y). 


Then  from  (3.9) 


a‘W+(l-12m)W=~^9"U(t'V)' 


Thus 


Sup(t,y)  «  ^dyu(t,  y)  -  f  e  {t  s),Td yu(s,y)ds 


-  -  I  e 
T  .) 


-{t-s)/T 


1  l 

2  (dvu(t>  v)  ~  dvu(s ,  y))  +  -  (d yu(t,  y)  -  d yu(s,  y))2 


-dyu(s,  y)ds. 


From  this  last  equation  we  are  persuaded  to  make  the  approximation 
Sup(t,y)  =  1 dyu(t,y )  -  ^  f  e~(t~s)/Tdyu(s,y)ds. 


(3.10) 
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4  The  dynamics  equation 

We  proceed  to  write  the  dynamic  equation  for  a  particular  rubber  material  considered  in 
[11].  For  this  particular  material  the  energy  density  proposed  in  [11]  consists  of  two  parts, 
Wcc  and  Wpc,  where  Wcc  is  the  contribution  to  the  energy  density  of  a  unit  CC-box  from  the 
cross-linked  molecules  constituting  the  box.  A  contribution  to  the  energy  density  function 
from  entrapped  PC-molecules  is  denoted  Wpc.  The  Wcc  and  Wpc  proposed  in  [11]  are  given 
by 


Wcc  —  105  (A^c  +  X2c  +  X\c  —  3)  + 

103  (X21c\22c  +  X21cX23c  +  \22c\23c  -  3) 

Wpc  =  169  (A^p  +  X22p  +  A 3p  —  3)  +  0.0138  (Alp  +  A^p  +  X3p  —  3)  + 

7.89  {X2lpXlp  +  A^pA 3p  +  XlpXlp  —  3)  , 

where  the  specihc  numbers  in  these  expressions  are  derived  from  validation  with  experimental 
data  of  Young  and  Danik  (see  [10]).  Using  the  fact  Ajc  =  and  X2c  =  £2  and  (3.1),  we  find 

W^cc  =  208 Up  +  OMp  +  •  •  •  . 


Thus 


dWCl 

duv 


=  416-Up  +  a2u2  +  a3ul  H - =  ge(uy )  «  416u 


v 


(4.1) 


Next  using  (3.3)  we  have 


WPc  =  169  (  X2lp  +  -  3  \  +  0.0138  ^Ajp  +  -  3 


dWr 


pc 


dX 


1 V 


'lp 


7.89  (  2Ai p  +  — ^ - 3 

^1  p 


'ip 


1  p 


=  169  (  2Alp  —  — p  )  +  0.0276  (  Xfp  +  -  3  )  [2Xlp-^ 


lp  X 


1  v 


X 


1  v 


23.67  2Alp  +  -y  -  3  2  -  —  )  . 


'lp 


Recalling  Aip  ~  1  +  Sup  we  have 


K  IQHSup. 

d  X  ip 
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Using  (3.5)  we  have 


dWpc 

dWpC  d\\p 

dXic 

d\\p  dXic 

Aip 

1  ~\~  5 Up 

Aic 

X\ c 

dWpc 

dwpcd  Alc 

dUy 

dXic  duy 

=  512  5up  ( 1 

\ 

~  512 5up. 

11  1 

1  +  5up  -  -dyu  -  -SupdyU  +  -(<9 yu)2  +  . . . 

/  11  1 
~  1014 5up  (  1  +  $up  ~  - dyii  —  -SupdyU  +  ~(dyu)2  + 


1 

2 


2 SllpdyU  +  ^(dyll)~  +  .  .  .^ 


Thus  to  first  order  we  have 


^  «  512 5up. 

OIL, 


(4.2) 


Now  we  consider  the  shear  vibration  of  a  rubber  sample  with  a  square  cross  sectional  area 
Ac  and  density  p.  From  (4.2)  and  (3.10) 


dW. , 


pc 


256 


du. 


256(9 yu(t,y) - /  e  s^Tdyu(s,y)ds. 


T 


(4.3) 


The  shear  dynamics  of  the  rubber  sample  are  given  by 


pAcd2u  -  dy 


8Wrr  DW, 


du. 


+ 


pc 


du 


=  F, 


where  F  is  the  external  load.  Using  (4.3)  and  only  first  order  terms,  we  have  our  desired 
linear  equation 


256A,  rt 

Jo 


pAcd2u  —  672 Acd2u(t,y)  +  /  e  0  s^r(9 2u(s,y)ds  =  F. 


T 


(4.4) 


5  A  nonlinear  reptation  model 

In  the  previous  section  the  deformations  of  the  PC-molecules  and  the  CC-moleculcs  were 
assumed  to  be  related  according  to  the  formula 

Ad  _  l 
AL  _  L 


leading  to  the  relationship 


f)  \  \ 

u/xjp  _  1 TPr 

d\ A,-  13 


He  'He 

between  the  stretches. 

The  above  relationship  has  the  merit  that  it  would  make  it  feasible  for  the  PC-molecule  to 
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remain  within  the  constraining  tube  during  relaxation.  However,  a  more  realistic  and  general 
relationship  is  of  the  form 


Ai 

A L 


f(t ,  j),  or  =  /(^ 

J—j  AjV* 


(5.1) 


where  /  is  a  material  dependent  nonlinearity.  With  this  more  general  relationship  between 
the  deformations,  the  relaxation  process  after  an  instantaneous  step-strain  in  the  time  inter¬ 
val  0  =  to  <  t  <  t\  is  given  by 


£(t)  —  fu  +  f(to,  —)AL0e  t 


In  the  time  interval  >  1  the  relaxation  process  has  the  form 


m—  1 


t(t)  =  Ip  +  y:  f(ti,  j-)A Lie 


ti) 


2—0 


In  the  limit,  as  At  —  tt  —  tends  to  zero,  we  have 


.  /‘f  Hs)  dL(s)  (*-») 


In  differential  form  we  have 


..  £(t)  dL  1  .  . 


(5.2) 


We  can  write  (5.2)  in  terms  of  the  stretches  of  the  PC-moleculcs  =  1,2,3  and  those 
of  the  CC-molecules  A ic,  i  =  1,2,3.  When  writing  (5.2)  in  terms  of  the  stretches,  we  take 
£0  =  1.  In  particular,  setting  e1  =  Sup  =  \\p  —  1,  we  have  from  (5.2)  the  equation 


£i  H — £i  —  Aic/(f,  — - - ). 

Aic 

In  the  problem  at  hand  A^A^  =  1,  and  from  (3.2), 

=  l-£  +  £2-^£3±---, 


(5.3) 


where  we  set  e  =  -u(y).  Next  we  write  (5.3)  as 

1  3 

£\  -( - £\  =  £{\  +  2e  +  — £2  +  •  1 

r  2 


■/(t,  1  +  £i  —  £  —  ,£i£  +  £2  +  £i£2  —  —  £3  A  •  •  • 


In  the  case  /  is  independent  of  t,  which  we  assume  henceforth,  we  expand  /  about  the  point 
1  and  obtain 


1  3 

£\  H - £\  =  d(l  +  2e  +  —  £2  +  •  •  • 

r  2 


1 


/(l)  +  /( l)(£i  ~  ~  £\£  +  £2  +  £i£2  —  -£3  + 


In  Section  3,  f(t,  — 1 — )  =  — 1 — .  In  this  case  (5.4)  becomes 

Ale  Ale 

1  3  n  , 

£\  H — £\  —  d(l  +  2e  +  —  £  +  •  •  • ) 
r  2 


1  ~h  —  E  —  E\E  -b  £  +  E\E  —  ~E  it  •  •  • 


We  note  that  a  special  case  of  this  latter  equation  is 

1  3 

£\  H — £\  —  d(l  +  2e  +  —  e  ) 
r  2 

=  5<£  +  £2  +  k)’ 

which  is  itself  a  special  case  of  the  general  model 

£l  H - £1  =  —(bi£  +  b2£^  +  b3£3). 

t  at 

These,  of  course,  are  both  special  cases  of  the  general  model 


)  + 


(5.4) 


(5.5) 


(5.6) 


£i  H — £i  —  9v(£,  £)■ 

T 

Neglecting  all  higher  order  terms  in  (5.5),  we  find 


(5.7) 


£\  H - £i  —  £, 

T 

which  is  (3.10)  of  Section  4. 

The  dynamics  of  the  rubber  (see  Section  4)  are  given  by 

,  o2  d  ,d\Vcc  ,  dwpc^  ^ 
pAcdt  u  —  w-(yy—  +  w— )  =  F. 


Now 


dy  duy  duy 


dW, 


pc 


dWpc  dXlc 


du 


dXic  duy 
dWpc  dXip  dXlc 
dXVp  dXXc  duy 
dWpc  , Aip.  1 
<9 A,p  Alc  2 

512£i/(1  T  £\  —  £  —  £\£  T  £  3”  £i£  —  —  £3  i 


(5.8) 


(5.9) 
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From  (5.9)  we  have  now 


pAcd2u  —  dy{ge(e )  +  512ei/(l  +  £\  —  e  —  E\£  +  £2  +  £\£2  —  -£3  ±  ■■■)}  =  F,  (5.10) 

where  ge(s)  =  832e  +  a-2£ 2  +  •  •  •  .  This  is  a  general  nonlinear  model  based  on  the  nonlinear 
constitutive  relationship  involving  /  of  (5.1)  and  that  of  (4.1).  In  Section  3  we  had  the 

simpler  proportional  assumption  /(— ^)  =  — ^  as  well  as  ge{£)  ~  832e.  In  that  case  (5.10) 

Ale  Ale 

becomes 

pAcd2u  —  (9y(832e  +  512ei(l  +  £i  —  £  —  E\E  +  £2  +  £\£2  —  H - ))  =  F.  (5.11) 


If  we  retain  only  first  order  terms  in  e,£\  in  (5.11)  we  obtain  the  equation 


pAcd2u  -  dy( 832e  +  512ei)  =  F. 


Thus,  a  linearization  of  (5.11)  and  (5.5)  gives 

pAcd2u-dy( 832£  +  512ei)  =  F 

1 

£\  H - £i  —  £, 

T 


which,  recalling  (3.10)  and  £ 


we  see  is  the  same  as  (4.4)  of  Section  4. 


6  The  nonlinear  model  and  simple  shear  experiments 

Model  development  and  identification  processes  (inverse  problems  with  experimental  data) 
for  nonlinearities  and  parameters  were  carried  out  for  filled  elastomers  undergoing  simple 
shear  deformation.  A  series  of  dynamic  experiments  were  designed  and  carried  out,  each 
involving  one  or  more  “double  sandwich”  fixtures  with  layers  of  filled  rubber  at  the  interfaces 
as  depicted  in  Figure  2. 


Load  Cell 


Figure  2:  Schematic  of  the  double  lap  simple  shear  experimental  device. 
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The  side  bars  were  fixed,  while  the  middle  bar  was  either  perturbed  by  an  impulsive 
hammer  hit  or  provided  with  an  initial  strain  and  then  released.  In  all  experiments,  both 
accelerometer  and  load  cell  data  were  collected  for  the  motion  of  the  middle  bar.  The  samples 
consisted  of  A225  highly  filled  rubber  with  dimensions  0.1  x  0.1  x  0.8  in.  with  two  such  layers 
in  each  of  the  double  lap  devices  as  seen  in  Figure  2. 

For  simple  shear,  the  model  developed  and  used  is  based  on  (5.10)  with  Kelvin- Voigt 
damping  added  as  given  by 


d 

pAcUtt  V '  D^tyy  7^  AC(J  0, 


0  <y< 


(6.1) 


{Muu(t,£)  +  CDuty(t,£)  +  Aca(t,£)}  =  F(t)  +  Mg,  (6.2) 

u(t,  0)  =  0,  (6.3) 

u(0,y)  =  %  ut(0,j/)  =  0.  (6.4) 

Here  p  is  the  mass  density,  Ac  is  the  cross  sectional  area  of  the  sample,  Cd  is  a  damping 
coefficient,  M  is  the  mass  of  a  tip  body  possibly  attached  to  the  middle  bar  (M  =  0  if  the 
tip  mass  is  not  included  in  the  experimental  set-up),  F(t)  is  an  applied  (at  y  —  i)  external 
force  when  present,  and  a  is  the  elastic  and  viscoelastic  shear  stress  given  in  general  form  by 


cr(e,ei)  =  ge(e)  +  e1 

(6.5) 

£i  +  -£i  =gv(e,i), 

(6.6) 

as  developed  in  the  previous  sections  of  this  paper. 

For  each  experimental  data  set,  a  least-squares  minimization  problem  was  carried  out 
to  identify  the  unknown  parameters  q  in  the  nonlinearities  ge  and  gv  using  the  accelerom¬ 
eter  data  {a*}  (the  load  cell  data  was  not  used  in  the  optimization  procedure)  with  the 
corresponding  model  solutions  Utt(U,£;q )  in 

1  " 

j(q)  =  ^^2\ai~Utt(U,£;q)\2-  (6.7) 

i=  1 

The  resulting  optimized  model  was  then  used  to  generate  load  cell  values  which  were  then 
compared  to  the  corresponding  experimental  load  cell  data. 

Initial  experiments  involved  one  double-lap  device  with  an  attached  tip  mass  (to  produce 
deformations  in  shear  sufficient  for  observation)  which  was  excited  with  a  hammer  hit  on  the 
middle  bar.  The  experiment  was  performed  with  highly  filled  rubber  samples  with  different 
amounts  of  attached  extra  mass.  In  general,  we  found  that  the  data  collected  with  no  extra 
mass  on  the  sample  could  be  approximated  well  with  the  above  type  models  even  when 
the  viscoelastic  response  function  gv  is  assumed  to  be  linear.  This  is  not  surprising,  since 
the  maximum  strain  levels  in  these  type  of  experiments  were  below  10%.  However,  when 
extra  mass  was  involved  to  achieve  larger  deformations  and  strain  levels  (this  required  huge 
masses!),  the  experiment  did  not  provide  suitable  data,  since  additional  “tilting”  modes 
became  excited  and  these  could  not  be  accounted  for  by  the  above  one-dimensional  model. 
We  subsequently  redesigned  the  experimental  device  to  consist  of  a  hub  with  four  identical 
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arms,  each  containing  a  double  lap  device  as  pictured  in  Figure  2.  The  arms  could  be  latched 
down  so  as  to  produce  initial  strains  of  up  to  100%.  Upon  quick  release,  this  device  produced 
very  clean  free-release  data  for  the  samples  undergoing  simple  shear  with  no  additional 
deformations;  these  proved  adequate  for  our  model  estimation  and  validation.  Thus,  the 
experiments  consisted  of  latching  down  the  four  arms  of  the  fixture  at  some  prescribed  initial 
strain  (50%,  70%  and  100%  strain)  and  then  suddenly  releasing  the  latch.  We  began  with 
approximately  50%  initial  shear  strain  and  repeated  each  release  experiment  three  times. 
We  again  collected  both  accelerometer  and  load  cell  data.  The  sequence  of  experiments  was 
repeated  with  about  70%  and  100%  initial  shear  strains,  respectively. 

For  our  model  fitting  we  used  a  number  of  different  constitutive  relationships  between 
the  stress  and  the  strain,  including  linear  and  nonlinear  functions  for  both  ge  and  gv.  The 
identification  problem  was  first  performed  for  the  data  obtained  for  100%  initial  strain. 
We  quickly  found  that  a  linear  viscoelastic  response  function  gv  was  no  longer  adequate  to 
describe  the  data  in  this  higher  strain  regime.  Thus,  we  were  led  to  models  of  the  form 


1 

£i  H — £i  —  gv 

T 


du  d  2u 
dy  ’  dtdy 


_  _  _  d£  1  d2u 

where  ge  and  gv  are  cubic  polynomials  and  gv  may  or  may  not  depend  on  —  =  -  -  , 

(If  2  dtdy 


re., 


gVi  =  gvd  if  the  nonlinearities  are  the  same  in  both  increasing  and  decreasing  deformation). 
The  notion  that  one  might  need  gv  ' s  of  the  form 


~  /  -x  j  9vi(uy)  if  iiy  >  0,  (■ uy  increasing) 

uy)  =  |  if  £  <  0_  (a'  decreasing),  <6-8) 

was  motivated  by  earlier  results  [4,  6,  7]  with  similar  materials  undergoing  tensile  deforma¬ 
tions  where  such  models  were  necessary  because  stress-strain  curves  clearly  indicated  that 
the  nonlinearities  are  distinctly  different  for  contraction  versus  elongation. 

The  parameters  in  a  linear  gvl  =  gvd  that  provided  the  best  fit  for  the  100%  initial  strain 

data  (see  Figure  3)  did  not  produce  a  model  that  described  well  the  70%  and  50%  initial  strain 

experiments.  Thus  we  next  used  the  full  nonlinear  hysteretic  model  (i.e.,  gvi  %  gV(i  in  (6.8) 
taken  as  cubic  nonlinearities)  to  approximate  the  100%  initial  strain  data  and  found  several 
different  sets  of  cubic  nonlinearities  that  gave  a  very  good  approximation.  Again,  earlier 
efforts  with  these  materials  suggested  that  cubic  (as  opposed  to  higher  order)  nonlinearities 
might  be  sufficient  to  capture  the  dynamics.  To  find  a  set  that  also  provides  a  good  fit  for 
the  70%  initial  data  we  optimized  on  both  data  sets  simultaneously  (i.e.,  we  used  both  data 
sets  in  the  cost  criterion  (6.7)).  The  best  fit  in  this  case  is  depicted  in  Figures  4  and  5.  We 
found  that  the  set  of  parameters  identified  in  this  case  described  both  data  sets  with  very 
good  accuracy  (Figures  4,5).  Finally,  we  found  that  the  50%  initial  strain  data  sets  were  the 
most  difficult  to  use  with  the  models  (see  Figure  6).  Our  investigations  suggest  that  in  this 
data,  the  deformations  are  likely  too  small  to  obtain  reliable  fits  due  to  unfavorable  signal 
to  noise  ratios. 
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Acceleration  data  and  the  approximated  values 


Comparison  of  load  cell  data  and  the  model  predicted  values 


Figure  3:  Hysteretic  model  with  gm  =  gvci  linear  for  100%  initial  shear  strain. 


Acceleration  data  and  the  approximated  values 


Figure  4:  Model  with  gvl  %  gv^  cubic  best  fit  compared  with  100%  initial  shear  strain  data 
response. 


13 


Acceleration  data  and  the  approximated  values 


Comparison  of  load  cell  data  and  the  model  predicted  values 


Figure  5:  Model  with  the  same  (as  in  Figure  4)  gvi  %  gvd  cubic  best  fit  compared  with  70% 
initial  shear  strain  data  response. 


Acceleration  data  and  the  approximated  values 


Figure  6:  Best  model  with  gvi  %  gVfj  cubic  compared  with  50%  initial  shear  strain  data 
response. 
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7  Analysis  of  the  linear  model  equation 


In  this  section  we  analyze  the  linear  model  equation  (4.4)  and  make  some  conclusion  about 
the  model.  Dividing  by  pAc  in  (4.4)  we  obtain  the  equation 


672 


256 


d?u-—d2u+—  /  e-<-t-'VT&u(s,y)ds  = 


P 


PT  Jo 


pA, 


F. 


We  will  later  on  solve  (7.1)  with  boundary  and  initial  conditions  given  by 

du 

u(t,  0)  =  0,  — ( t,b)=G(t ),  u(0,y)  =  dtu(0,y)  =  0. 


To  facilitate  our  analysis,  we  set 


672  256  I'1 

v  = - u - 

P  PT  Jo 


D-(t-s)/r 


u(s,  y)ds. 


Then  from  (7.1)  we  obtain 


256 


416 


vtt 


-vt 


672 r  672 r2 


672 

- 1 

P 


Jyy 


re-*™ 


672 t 
672 


vs  H — v  I  ds 


T 


672  (  416 


p2Ac  p  V672 r. 

We  rewrite  the  initial  and  boundary  conditions  in  (7.2),  using  (7.3),  as  follows: 


(7.1) 


(7.2) 


(7.3) 


F - e-^‘«(0,y).  (7.4) 


P 

464 


^  672r  ^ 


v(t,  0)  =  0 

vSy(s,b)  +  -Vy(s,b) 
T 

v(0,y)  =  o 


ds  +  dyu(0,  b)e 


vt(0,y)  +  ~v(0,y)  =  0. 

T 


G(t) 


(7.5) 


The  second  equation  in  (7.5)  can  be  written  as 


P 

672 


vy(t,b ) 


—  41(3  +  ,  \ 
-e  672rr^(0,6) 


=  <?((). 


(7.6) 


We  next  consider  the  Sturm-Liouville  problem 


-Z”~\Z  =  0 
Z(a)  —  h0Z\a)  =  0,  h0  >  0 
Z(b)  -  luZ'ib)  =  0,  hi  >  0. 


We  recall  that  there  is  a  corresponding  sequence  of  eigenvalues  0  <  Ai  <  A2  <  . . .  <  An  < 
...  y  00  and  corresponding  orthonormal  eigenfunctions  (pi,  </?2,  •  •  • ,  <pn, _ 
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Consider  a  solution  of  (7.4)-(7.6)  in  the  form 


OO 

v(t,y)  =  ^Vn(t)(Pn(y)- 

71—1 


Writing 


where 


f{t,y) 

™lf 

P2AC 


OO 

^fn(t)(Pn(y), 

71—1 


672  /  416 
~p~  \672r 


2 

416  f  .  . 

e  672rC^(0,l/), 


we  have  from  (7.4)  the  equation 

.  256  , 

-„(*)  + 

2  pt 

Jo 


( 672  ^ 
+  — A 
V  P 


416 


A 


416 

672r 


416 
g  " 672t 


672r2 )  n 

t‘~‘>  (hM  +  h(»)l  *  =  /»(*)■ 


(7.7) 


Set 


wn{t) 


Then 


w'n{t) 


Let 


From  (7.7)  we  have 


d 

dt 


416  \ 
6727/ 


2  pt 

J  e'^M  ^w(s)  +  -wn(s 


fis. 


416  \  /  ,  ,  ,  1  .  A  416  ,  . 

era?  (”»W  +  -”»W 


<(()  = 
<(/)  = 

<(/)  = 

d  = 

An  = 


^n(0 

v'jt) 

wn(t) 

416 

672r 

672 

- A?. 

P 


416 
672r2 ' 


/  0  1  0  \ 

A  -256 

-An  -  -1 

/  672t 

-  d2 

T 


-d 


(7.8) 


(7.9) 
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Next  define 


An  — 


(  0 

-A, 


1 

-256 


rf.  672T 

-  d 2 

T 


0  \ 

-1 


-d 


Then 


\^i-An\  =  e+ 


(  256 


V 


672r 


+  d  £2  +  (  d2  + 


256 


d 2 


672rrf  +  An  )  ^  +  +  V' 


The  Routh-Hurwitz  Theorem  guarantees  that  the  roots  of  the  polynomial  equation 

z ^  T  o,\z  T  (X2Z  T  03  —  0 

have  negative  real  parts  if  the  following  conditions  are  met: 

(i)  >  0 


(ii) 


CL\  CI3 

1  a2 


>  0 


(iii)  <23  >  0. 

Thus  the  roots  of  | £/  —  An  |  =  0  have  negative  real  parts  if 


256 

672r 

1 


d 


.  d2 
And  T  — 


256 

672r 


r 


d  T  d?  T  Ar: 


>  0. 


Thus,  the  eigenvalues  of  7Ln  have  negative  real  parts  if 

2 


d 


f 


256 


V  672r 


+  d  )  + 


/  256  \  .  d2 

\672r  +  d)  An  -  -  ~  >  °- 


(7.10) 


Now  recalling  that 


d  = 
A  „  = 


416 

672r 

672 

/\j 

P 


416 
672r2 ’ 


a  simple  calculation  reveals  that  the  inequality  in  (7.10)  is  true  if  Xn  >  0.  Thus  the  eigen¬ 
values  of  An  have  negative  real  parts. 

From  (7.8),  (7.9)  we  have 


vn(t) 

v'n(t)  ]  =  e“" 
Wn{t) 


Vn{ 0)  \  rt  /  0 

%»(0)  +  /  e^s)An  fn(s)  I  ds. 


0 


0 


Remark  7.1:  Since  the  eigenvalues  of  An  have  real  parts  negative  we  have 

fn(0) 

«n(0) 

0 


. .  t  .4 ,1 


as  t 


00. 
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8  Solution  of  the  linearized  model  equation 


In  this  section  we  proceed  to  solve  (7.4),  (7.5).  For  convenience  we  rewrite  (7.4)  and  (7.5) 
as  follows: 


256  416 

Vtt  +  672rVt  ~  672 r2 


v  — 


672 

P 


-yy 


(  416 
\672r 


672 

P2K 


'  o 

672 


- 416_ / ^ — \ 

^  672r  v6  *7 


416 


vs  H — v  )  ds 

T 


e  (8-!) 


The  initial  and  boundary  conditions  are 


v(0,y)  =  0, 


vt(0,  y)  =  0,  (8.2) 

v(t,  0)  =  0, 


and 


P 

672 


vy(t,b ) 


—  4it)  +  /  „  7  x 
—  e  672r^  vy(0,  b) 


+ 


256p  f 


(672)2r  J0 


e  672 r(f  S)ny(s,6)ds 


+G(0+)e-6W4  =  G(t). 


(8.3) 


Set 


w  = 


672r^2/(0,5)  + 


256p 

(672)2r 


416 

672r 


^  S\>y(s,b)ds. 


(8.4) 


Then  equation  (8.3)  becomes 


g^y(t,  6)  +  w  +  G(0+)e  Aw®*  =  G(f). 
Differentiating  (8.3)  we  have 


(8.5) 


p  d 
672  df 


fyM) 


416 

672r 


w  + 


256p 

(672)2r 


416 

672r 


G(0+)e-^i  =  G’(t), 


(8.6) 


where  we  have  assumed  that  G  is  sufficiently  regular,  e.g.,  G  G  C1.  From  (8.5)  and  (8.6)  we 
have 


p  d 
672  dt 


Vy(t,b)  + 


P  vy(t,b)  =  G'(t)  +  2^-_G(t), 


672 r 


672r 


(8.7) 


and  from  (8.2)  we  see  that 

vy(0,b)  =  0.  (8.8) 

From  (8.7),  and  (8.8)  we  thus  see  that  we  can  give  vy(t,b )  explicitly  in  terms  of  the 
function  G.  Thus  we  write 

vy(t,b)  =  rj(t) 
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where  rj  can  be  explicitly  determined  in  terms  of  the  function  G.  Thus,  we  replace  (8.2)  and 
(8.3)  by 


v(0,y)  =  vt(  0,y)=0 

v(t,  0)  =  0  (8.9) 

Vy(t,b)  =  rj(t). 


We  now  proceed  to  solve  (8.1)  with  initial  and  boundary  conditions  given  by  (8.9).  To¬ 
ward  this  goal  let 

Q(t,y)  =  v(t,y)  -yr](t). 

In  terms  of  Q  the  equation  (8.1)  becomes 


Q 


t.t 


256 

672r 


Q,  -  ««  Q  - 
V  672r2 V  p  V 


yy 


416  \2 
672r  ) 


e-e^h-s)  lQs  +  tQ)ds  =  f,  (8.10) 


where 

~  „  256  .  416  .  .  /  416  \  ft  416  r-L  \  f  /  1  1 

f  =  f-yv  (t)  -  ^vn  ( t )  +  M2^yv(t)  -  (^— J  ^  e  j  j rn  (s)  +  -yv(s)jds. 

For  initial  and  boundary  conditions  we  have 

<5(0, 2/)  =  —2/77(0) 

Qt(o,y)  =  -W(  0)  (8.11) 

Q(t,  0)  =  Qy(t,  b)  —  0. 

Let  us  now  consider  the  Sturm-Liouville  problem 

~z"-\z  =  0 

Z{  0)  =  0  (8.12) 

Z\b)  =  0. 


The  eigenvalue  problem  (8.12)  has  a  system  of  eigenvalues  /b  <  /b  <  ■  ■  ■  S'  00  and  a 
corresponding  orthonormal  system  of  eigenfunctions  tp]_,  ^2)  ■  ■  ■ ,  V’nj  ■  ■  ■■  Let 


R(t,y)  = 


/  416  \ 
\672r  ) 

J'  e-^{t~s)  Iqs 

R(t,y ) 

00 

=  ^2rn(t)ipn(y)t 

n=  1 

f(t,y) 

00 

=  ^fnWniy), 

71—1 

Q{t,y) 

OO 

=  ^2qn(t)lf>n(y)t 

71—1 

19 

1 

-( 

r 


We  now  write 


and  seek  a  solution  for  (8.10)  and  (8.11).  Let 


where 


Then 


where 
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672r 

d2 
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d2 

— d 
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672  416 

~V^n  ~  672^’ 


(  qn(t) 
\  rn{t) 


(  9n(0)  \  rt  (  0 

=  (  q'n(0)  I  +  j(  f  fn(s ) 


ds, 


Qn(  0) 
9n(°) 


-V(Q)  /  y^n(y)dy 


-V(0)  /  y^n(y)dy. 


Once  Q  is  determined  the  solution  v  for  (7.4)  and  (7.5)  is  given  by 


v(t,y)  =  Q(t,y)  +  yrj(t). 

Remark  8.1  The  eigenvalues  of  An  have  negative  real  parts  if  yn  >  0.  This  is  established 
using  the  Routh-Hurwitz  theorem  as  was  done  earlier. 

Remark  8.2  It  can  be  shown  that  the  eigenvalues  of  An  are  such  that  their  real  parts  are  less 
than  a  fixed  negative  real  number  for  all  n. 

Remark  8.3  Because  all  the  eigenvalues  of  An  have  real  parts  negative,  we  have 


f  ?n(0)  \  /  0  \ 

(  9^(0)  I  -»■  (  0  I  ,  as  t  ->  oo. 


9  Concluding  remarks 

The  stress  strain  models  of  (6.5)-(6.6)  given  by 


cr(e,£i)  =  ge{e)  +  e1 

1  /  *\ 

El  H £  1  —  gv{£,  £), 

T 

discussed  in  Section  6  are  special  cases  of  more  general  models  in  which  one  admits  hetero¬ 
geneity  in  the  PC-molecules.  Specifically,  one  might  have  a  distribution  of  different  classes 
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of  PC-molecules,  each  characterized  by  different  internal  strains  £in  depending  on  physical 
properties  such  as  length,  relaxation  time,  etc.  For  example,  the  PC-molecules  might  have 
differing  relaxation  times  r  so  that  =  £in(t,y ;  r).  For  a  continuum  one  would  then  have 
a  probability  distribution  P(r)  on  r  G  T,  where  T  is  the  set  of  possible  relaxation  times. 
The  corresponding  stress  strain  models  are  then  given  by 

V\  P )  =  9e(e(t,  y),e(t,  y))+'r  ein(t,  y\ r)dP(r),  (9.1) 

where  £in(t,y ;  r)  satisfies,  for  each  r  e  T, 

£in(t,y;r)  +  -ein(t,j/;T)  =  &,(e(t,y),e(t,j/)).  (9.2) 

r 

The  models  we  used  with  shear  data  in  Section  6  required  only  a  simple  Dirac  delta  measure 
dP(r )  =  8T1(r)dT  with  =  eq  then  satisfying 

H - £i  =  fik(£>  P) 

T\ 

in  order  to  obtain  reasonable  fits  to  the  experimental  data.  We  compare  this  with  the  situ¬ 
ation  for  tensile  deformations  reported  in  [4,  5]  where  the  data  required  multiple  relaxation 
times  in  the  model  so  that  the  measure  P  had  atoms  at  7q  and  t2,  (i.e.,  the  measure  is 
composed  of  Dirac  measures  concentrated  at  iq  and  t2),  and  the  constitutive  law  led  to  the 
model 


<r{t,V,P)  =  9e(£(t,y),£(t,y))+-/1£1(t,y,T1)+~/2£2(t,y]T2) 

1 

£\  +  — £\  — 

Ti 

1 

£2  +  —^2  = 
r2 

In  the  shear  studies,  we  found  that  the  nonlinear  models  (specifically,  gv)  were  essential  in 
describing  our  data,  but  a  simple  uniform  relaxation  time  was  adequate  for  the  PC-molecules. 
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